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The magnetic field dependence of energy levels in gapped single- and bilayer graphene quantum 
dots (QDs) defined by electrostatic gates is studied analytically in terms of the Dirac equation. Due 
to the absence of sharp edges in these types of QDs, the valley degree of freedom is a good quantum 
number. We show that its degeneracy is efficiently and controUably broken by a magnetic field 
applied perpendicular to the graphene plane. This opens up a feasible route to create well-defined 
and well controlled spin- and valley-qubits in graphene QDs. We also point out the similarities and 
differences in the spectrum between single- and bilayer graphene quantum dots. Striking in the case 
of bilayer graphene is the anomalous bulk Landau level (LL) that crosses the gap which results in 
crossings of QD states with this bulk LL at large magnetic fields in stark contrast to the single-layer 
case where this LL is absent. The tunability of the gap in the bilayer case allows us to observe 
different regimes of level spacings directly related to the formation of a pronounced "Mexican hat" 
in the bulk bandstructure. We discuss the applicability of such QDs to control and measure the 
valley isospin and their potential use for hosting and controlling spin qubits. 

PACS numbers: 73.21. La, 81.05.Uw, 74.78. Na, 71.70.Di 



I. INTRODUCTION 

Graphene is one of the most promising materials for fu- 
ture nano-electronics.^'^ This is related to its truly two- 
dimensional character yielding perfect electron confine- 
ment in one spatial dimension. In order to build func- 
tional nano-devices such as single-electron transistors, 
quantum point contacts, and quantum dots (QDs), ad- 
ditional confinement in the remaining two spatial dimen- 
sions is needed. Due to the absence of a gap in the spec- 
trum, this is a rather demanding task in both single- and 
bilayer graphene, in contrast to electrostatically defined 
QDs in semiconductors such as GaAs. One possibility of 
overcoming this difficulty consists in etching or scratching 
nanostructures into graphene flakes. This has been done 
to experimentally study, for instance, transport through 
graphene nanoribbons,'^'^'^ single-electron transistors,^''' 
and, very recently, even QDs showing pronounced sig- 
natures of excited states.^ Nevertheless, to increase the 
functionality of graphene nano-devices it is desirable to 
develop gate-tunable structures. 

In this article, we study the energy spectrum of gatc- 
tunable QDs both in single-layer and bilayer graphene. 
In single-layer graphene, we assume a constant gap in the 
whole system that might be introduced by the underly- 
ing substrate. ^'^"^ In bilayer graphene, it is well-known 
that a gap can be generated by applying different elec- 
trostatic potentials to the upper and lower layer, ^^'^^ 
which has already been experimentally observed. ^^'^^'^^ 



Once there is a physical mechanism that gives rise to 
a gap, bound states exist in the presence of an electro- 
static confinement potential. We focus on the magnetic 
field dependence of bound states in circularly symmet- 
ric QDs. Whereas previous work has analyzed bound 
states in single-layer graphene subjected to spatially in- 
homogeneous magnetic fields, we analytically study the 
magnetic-field dependence of bound states due to elec- 
trostatic (i.e. non-magnetic) confinement. A comple- 
mentary numerical analysis has been done to study the 
Fock-Darwin spectrum of parabolic QDs in single-layer 
graphene^^ where only quasi-bound states but not true 
bound states exist. ^® 

Most remarkably, we show how the valley degeneracy 
can be lifted by an external magnetic field applied per- 
pendicular to the surface. This is of particular impor- 
tance to form valley-filters, -valves, or -qubits, and 
spin qubits^^ in graphene. To do so, it is essential to 
have full control both over spin and valley degrees of 
freedom and we show that a magnetic field is all that 
is needed to achieve this goal. Some of us have demon- 
strated that such a control can also be achieved in single- 
layer graphene ring structures with an Aharonov-Bohm 
fiux applied. Here, the emphasis is on the more feasible 
situation of a constant magnetic field applied to the whole 
system. We would also like to mention that the bro- 
ken valley degeneracy has an interpretation in terms of a 
magnetic moment that depends on the valley isospin. 
Our models for single- and bilayer graphene QDs are ap- 
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propriate for the physical situation of a smooth crossover 
between the dot region and the barrier. Therefore, atom- 
ically sharp edges do not play any role in our analysis 
which seems to be the most relevant case for possible ex- 
perimental realizations of gate-tunable QDs in graphene. 
Indeed, in Ref. 7 the absence of a four-fold level degener- 
acy in the transport data — due to valley and spin degrees 
of freedom — was attributed to intcr-valley scattering at 
the atomically sharp edges whereas the absence of spin 
degeneracy could result from spin scattering at dangling 
bonds at the edge of the QD. Both possible sources for 
an uncontrolled lifting of degeneracies are not relevant for 
our QD realizations. 

In the regime of strong magnetic fields, the bound 
states of single- and bilaycr graphene merge into the ap- 
propriate bulk Landau levels (LLs) as expected. How- 
ever, the nature of these LLs are quite different for the 
cases of single layer graphene versus bilayer graphene. In 
the bilayer case, one of the LLs crosses the gap with in- 
creasing magnetic field. QD bound states can cross this 
LL at large magnetic fields which leads to certain con- 
straints to form operational QDs. 

The paper is organized as follows. In Sec. II, we dis- 
cuss our model and the results for single-layer QDs. In 
Sec. Ill, we treat the bilayer case, in Sec. IV we discuss 
possible applications of our results for the emerging field 
of valleytronics and to spin-based qubits in graphene, 
and, in Sec. V, we draw our conclusions. 

II. QD IN SINGLE LAYER GRAPHENE 



Fig. 1). We also include a homogeneous magnetic field 
B perpendicular to the graphene plane. 

The Hamiltonian in the valley-isotropic form is given 

by23 

Hr = Ho+TAa, + U{x,y), (1) 

where Hq = v{p + eA) ■ cr, B ^ V x A = (0,0, B), 
V = 10^ m/s is the Fermi velocity and r = ± differenti- 
ates the two valleys K and K'. We choose the symmetric 
gauge A = Y(~yi2;,0) and assume a circular symme- 
try in the confinement potential U{x,y) = U(r) with 
r = y'x^ + y^. The vector operator cr acts on the A, B 
sublattice components of the spinor wave function and 
its vector components are given by the standard Pauli 
matrices. 

may be transformed into polar coordinates [(x, y) = 
(rcos(p,rsin(p)] (with h=l) 

Since Ht commutes with the total angular momentum 
operator J ^ = —id^ + azj^-, the energy eigenspinors can 
be chosen to be eigenstates of Jz 

^^(^,.)=e^™(,f;i.), (3) 

with i the eigenvalue of which has to be an half-odd 
integer. 
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FIG. 1: QD in single-layer graphene with a constant mass 
term A. An electrostatic potential with height Uo gives rise to 
bound states (dashed line) in the conduction band (c) defining 
a QD of radius R. Note that the confining potential U{r) is 
repulsive for holes in the valence band (v). 

In this section we study graphene in the presence of a 
constant mass term A (inducing a gap 2A) that might 
be introduced by the underlying substrate. The QD 
is defined by gates introducing an electrostatic confin- 
ing potential for electrons in the conduction band (see 



A. Bound state solutions 

To solve the eigenvalue problem Hr'^'^{r,(p) = 
E'i>'^ (r. If) we have to analyze 



with x^(f) = {XAi'^),XBi'^)V and 

Ht{i") = —ivdxdr + tAcTz + U{r) + 

+ 



(4) 



va-y 











3 + 1/2 



eBr 



(5) 



First, we solve Eq. (4) with a constant U{r) = Uq. Defin- 
ing e = E — Uq and b = eB/2, we obtain the following 
decoupled second order differential equations 

r^dhUr) + rdrXlir) - {b'r* + + nl)xl{r), (6) 

with (T = ±1, the upper sign corresponding to the A, the 
lower to the B sublattice, the coefficients entering Eq. (6) 
can be expressed as 



n„ 



2h{j + a/2) - (e^ - A')/v\ (7) 
|j-^/2|. (8) 
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Note that Eq. (6) does not depend on the valley index r 
anymore. However, Xai''^) depends on r through Eq. (4). 
The solutions to Eq. (6) are the confluent hypergeomet- 
ric functions M{a,b,z) and U{a,b,z). The boundstate 
solutions for the QD have the form 



ac,U{qa, 1 



,br^ 



r>R, 



P„M{q^,l + n„,br^) , r<R, 



(9) 



where = + 2(1 + n^)]. Eq. (9) is the general 

solution for waves that are regular at the origin and which 
decay exponentially as r ^ oo. We want to find the 
bound states for the following hard-wall potential 



C/(r) 



Uo , r> R, 
, r < R, 



(10) 



and define the corresponding energies as e< = E and 
e>=E- Uo. 

The ratios as/ctA and /Sb/Pa in Eq. (9) are fixed by 
the coupled first-order differential equation Eq. (4). This 
provides us with the general solutions for r < R and 
r > R. The matching conditions of the spinors at r = 
R gives then the eigenvalues and eigcnfunctions of the 
bound states. 



For j > 0, we obtain the following characteristic equa- 
tion for the allowed eigenenergies E of the QD 



^+ M (q< , J + 1/2, x) U{q> , J + 3/2, x) - ^+ M(g< , j + 3/2, x) U{q> , j + 1/2, x) = 0, 



(11) 



and for j < we obtain 

M(g<, -J + 3/2, x) U{q> - 1, + 1/2, x) - e< M{q< - 1, -J + 1/2, x) U{q>, -j + 3/2, x) = 0, (12) 



where x = bR'^={l/2){R/lBY with Is = ^h/eB the (Eq. (11)) becomes 
magnetic length. Without loss of generality, we choose B 
positive. The bound state levels for B negative can be ob- 
tained from the symmetry Hrij, B) = H^r{—i, ^B). We 
further introduced the parameters (/< > = (j — 1/2) 0{j) + 
1 - (e|,> - A2)/46i;2, i+ = (e< - TA)/4(i + 1/2), 

it = W(^> + rA), e< = U - l/2)/(e< + tA) and 

C> = l/(e> + tA) with 0{x) the Hcavisidc function. x Jj+1/2 ( ^ Je^ - A^ ) ifj_i/2 ( - 



In the limit of small magnetic fields (a; <C 1), the hyper- /2 , \ /2 

30metric fun 
13 in Rcf. 24) 



geometric functions reduce to Bessel functions (see Ch. + Jj-1/2 [ Jy ^< ~ ] ( 




(15) 



M(g<,n,x)=r(n) (-xq^f-''^/' 



X J„_i (2^y-xq<) , (13) and r 



For j < (Eq. (12)), we obtain Eq. (15) with j 



-T. 



Even in the limit of zero magnetic field, the character- 
istic equation (Eq (15)) cannot be solved in closed form 
in general. However, the fact that E{j,T) ^ E{—j,T), 
but E{j, t) — E{—j, — r) lies at the heart of our current 
2 approach to control the valley degeneracy by a magnetic 

J7(g>,n,a;) = (x(7>)'"^~"^^^ field. The first statement is a consequence of effective 

r(l + — n) time-reversal symmetry (eTRS) breaking within a single 

X A'„_i {2y/x g>) , (14) valley by a finite mass A.^s Formally, [Hr,f] ^ where 

T = iayC with C the operator of complex conjugation. 
The second statement is that the true TRS (which cou- 
where we have introduced the QD level spacing S = pies the two valleys) is not broken by a boundary alone, 
hv/R. For B = 0, the characteristic equation for j > i.e. at S = (see also subsection IV A). 
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FIG. 2: Bound state levels as a function of the QD radius R 
with (7o = A and for j = 1/2 at zero magnetic field. Full lines 
correspond to t = +1, dashed lines correspond to r = — 1. 



B. Results 

We first consider zero magnetic field. In Fig. 2 we 
show the energy levels of the QD as a function of the 
dot radius i?, evaluating Eq. (15) for j ~ 1/2. Full lines 
and dashed lines correspond to the two valleys. Due to 
the symmetry E{j, r) = E{^j, — t). the two set of curves 
display also the cases j = 1/2 and j ~ —1/2 in the same 
valley. The different solutions for the dashed and full 
lines are therefore a direct consequence of eTRS breaking 
in a single valley at zero magnetic field. However, if both 
signs of J were included, one would observe that the valley 
degeneracy was not broken at B = 0. 

In Fig. 3 we show the bound states of the QD as a 
function of magnetic field evaluating the characteristic 
equations Eqs. (11) and (12) numerically. In Fig. 3(a) 
we show the low-lying bound states in the conduction 
band. Note that the valley-degeneracy (or orbital de- 
generacy) is broken at finite magnetic field. The largest 
level spacing between the (non-degenerate) groundstate 
and first excited state we estimate from Fig. 3(a) to be 
at R/Ib ~ 1.8 and is about 165 mcV/i?[nm] for the pa- 
rameters used in Fig. 3. Ati?//B~1.8wc obtain for the 
Zeeman splitting = gf^sB ^ 200 meV/i?^[nm] using 
g = 2 which shows that the level spacing is always larger 
than the Zeeman energy for reasonable dot sizes. 

Considering a QD with i? = 25 nm, wc obtain a valley 
splitting Akj^i at R/Ib ~ 1-8 of about 6.6 mcV corre- 
sponding to 77 K, being much larger than 4 K, the tem- 
perature achieved by cooling with liquid helium. The 
necessary magnetic field corresponding to R/Ib = 1-8 
is B=3.41 T (and B = 0.85 T for R ^ 50 nm with 
Ak,k' ~ 3.3 meV) which is also easily achievable in 
the laboratory. A gap of size 0.23 eV has been con- 
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FIG. 3; a) Numerical evaluation of characteristic equations 
(11) and (12) as a function of R/Ib with Ib = {h/eBf^ the 
magnetic length and R the QD radius. We use A = 10 5 and 
f/o = A. a) The parameter regime of small B-fields where we 
observe a breaking of the level degeneracy. The full lines are 
for r = 1 and dashed lines are for r = — 1 corresponding to 
the two valleys of graphene. b) Same parameters as in a) , but 
for larger magnetic fields. The energy levels converge to the 
bulk Landau levels with increasing R/Ib- 



eluded from ARPES data in graphene on top of a SiC 
substrate. Therefore, the gap A and also the confining 
potential step height Uo could easily be larger than the 
QD level spacing S which is about 26 meV. These results 
suggest that such QDs confined in graphene would be an 
ideal host for spin qubits where the orbital degeneracy is 
controllable by a magnetic field. 

In Fig. 3(b) we show the merging of the QD states with 
the bulk Landau levels (LLs) 

En - ±S^{A/6y + 2niR/lB)^; n = 1, 2, 3, ... (16) 
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Note in particular, that 
— tA which Ues entirely 



with increasing magnetic field, 
there is a zero mode LL at ii^ = 
in one valley. 

In the next section we consider QDs in bilayer graphene 
where a voltage tunable mass gap is possible. 



III. QD IN BILAYER GRAPHENE 

In bilayer graphene, an electric field perpendicular to 
the layers generates a gap in the spectrum in a similar 
fashion to the staggered sublattice potential in the sin- 
gle layer. In this section, we will investigate the bilayer 
analogue of the single layer QD studied in the previous 
section, as shown in Fig. 4. Wc use the simplest nontrivial 



bilayer graphene 




dopant atoms 
./. . . 



SiOo 



back gate 



FIG. 4: QD in bilayer graphene: A back gate and dopants on 
top of the bilayer control the voltage V between the layers — 
leading to a controllable gap opening — as well as the Fermi 
energy (band filling) . An additional top gate allows to induce 
a spatially inhomogeneous electrostatic potential U (r) analo- 
gous to the single-layer model which leads to bound states in 
the conduction band (or valence band) of the bilayer. Another 
possibility is to use a split top gate (instead of a combination 
of top gate and dopants) to achieve a similar confinement. 

form of the Hamiltonian that captures the most impor- 
tant features of the spectrum and calculate the quantized 
energy levels of the QD as a function of the magnetic 
field and the relevant parameters of the band structure 
and the QD. The approximate Hamiltonian (we use) cor- 
rectly describes the crucial formation of an electronic gap 
in biased bilayer graphene. ^"^ We briefly discuss the issue 
of neglected terms in Sec. IIIC. 



A. Solving for the energy levels 

We separate the Hamiltonian in the bilayer into two 
parts: H = Ha + W[. Ho encodes the motion of the 
electrons within the planes and is given by two copies of 



the Dirac equation. In the valley-isotropic representation 
it takes on the form {h = v = 1) 



Ho — 



px + ipy 
Px - ipy 






Px - iPy 

Px + iPy 



(17) 



in both valleys. Like in the case of the single layer we 
add a magnetic field by the minimal coupling prescription 
P (p + eA) with A = {B/2){-y, x, 0). The other part 
of the Hamiltonian (i.e. HI) encodes the biasing field and 
the hopping t± between the two planes. The interplane 
hopping matrix element t±^ has recently been measured 
to be t± = 0.40 eV.^^'^^ In the simplest approximation 
we may take 



HI 



/ tV 
/ 2 



Vo 









tV 
2 






2 







TV I 

2 / 



+ C/(r)l, (18) 



with U{r) the applied electrostatic potential profile again 
given by Eq. (10). The index r = ±1 again distinguishes 
the two valleys (note that in the valley-isotropic repre- 
sentation the basis is chosen such that the two planes in 
the bilayer are exchanged in the spinors that describes 
different valleys). In Ref. 30, the same Hamiltonian H at 
zero magnetic field was used. However, the confinement 
described in Eq. (18) by U{r) was achieved in Ref. 30 by 
a position dependent "mass term" V{r), instead. 

To diagonalize H (i.e. to find the eigenspinors VE* that 
fulfill H'i> = E'i') we go to cylindrical coordinates in 
which the states are easily classified according to their 
conserved value of total angular momentum m (m being 
an integer). More explicitly, we factor out the angular 
dependence of the states according to 



/I 

e""'^ I e-'^P 

10 

Vo e"^y 



(19) 



Note that the angular momentum in the bilayer case is 
an integer m, in contrast to the half-odd integer j in the 
single layer case, which reflects the different pseudospins 
in the bilayer (pseudospin 1) and single-layer (pseudospin 
1/2). With the definitions j = m + 1/2 and s = sign(i?), 
the Hamiltonian Ho, which now acts on ^'i, can be writ- 
ten as 



Ho = 



I - (.7 - i)/e - si 

\ 



\ 



Se+j/e + ^e 
3ii ~si / 



(20) 
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In the latter equation, we have introduced the dimension- 
less coordinate ^ = r/{^/2lB), where Ib = ^h/{e\B\) is 
the magnetic length. The eigenvalue problem can now be 
solved by using the general solutions of the ordinary dif- 



J 



ferential equation imposed by T^o- The general solutions 
can be conveniently written in a simple way using the 
following functions (valid for all integers m and s — ±1): 



= e-«'/2^l"+"l+i/2^f([|m + a| + 1 + s(m - 1 - a)]/2 + k^/A, 1 + |m + aie)in^ + \m + a\). (21) 



These solutions are regular at the origin and are used 
for r < R. Note that k is an arbitrary parameter, 
which is chosen to be proportional to the energy eigen- 
value of the first subblock of the matrix in Eq. (20), i.e. 
Ho"^! = — iK^'i/V2/s for the first two components of 
(This choice is motivated by mathematical convencience 
to simplify the recursion relations in Eqs. (23a) — (23d) 
below.) In addition to Eq. (21), there are solutions that 
are irregular at the origin but vanish exponentially for 
r — > cx) which we use for r > R. These solutions are 
given by the same expression as Eq. (21) with the sub- 
stitution 

M{[\m + a\ + 1 + s{m - 1 - a)]/2+ 

kV4, l + \m + a|, e^)/r(l + Irn + a\) 

U{[\m+a\ + l+s{m-l-a)]/2+K^/A,l+\m+a\,^^). 

(22) 

For both types of solutions one can show the following 
identities by straightforward manipulations using the re- 
cursion relations for the confluent hypergeometric func- 
tions (see e.g chapter 13 of Ref. 24). 



and for r < R and m < — 1, 



(5e - {j - i)/e - 


= a{ ( 




(23a) 


{di + {j - l)/e + sO0m 


= a2< 


t'm-1, 


(23b) 


+j7e + sO0m+l 


= a3< 


Pm, 


(23c) 


(5« - j/i - sCl^t^m 


= a\ ( 


pm+1- 


(23d) 



For r < R and for m > 1 we obtain 

2 



01 = 

02 = 

ai = 



For r < R and m = 0, 



^72, 
2, 
2, 

(^2 -4s)/2. 



2 

^V2, 
2, 



(24a) 
(24b) 
(24c) 
(24d) 



(25a) 
(25b) 
(25c) 
(25d) 



ai 



2, 

«V2, 

(«' - 4s)/2, 
2. 



(26a) 
(26b) 
(26c) 
(26d) 



For r > R and all integer m we obtain 

al = -[{s+l) + n^{l-s)/i], 
4 = -[(l-.s) + Aj2(i + s)/4], 

al = -[(s+1) + (1 + kV4)(1-s)]. 
Therefore, by combining the solutions in the form 



(27a) 
(27b) 
(27c) 
(27d) 

















(pm-l 




























(pm+l 



*2 



(28) 



the Ti.0 part of the Hamiltonian (now acting on ^'2) can 
be replaced by: 



^^/2^i 



af 

; 

a§ 

. a| 0, 



(29) 



We now note that the transformations in Eq. (19) and 
Eq. (28) commute with the part of the Hamiltonian Til 
of Eq. (18). Therefore the task of finding the eigenvec- 
tors is transformed into the simple problem of finding the 
eigenvectors of a 4 x 4 matrix. Explicitly, the eigenvalue 
problem is equivalent to finding the non-trivial solutions 
of 
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— ia|/\/2/B 
t± 

\ 



E 



~ial/V2lB 
^ + U{r)-E 





t± 



^ + U{r)- 
~ia%l\f2lB 



tV 
2 



+ ;7(r) - £;/ 



^•2 = 0. 



(30) 



The non-trivial solutions are identified by finding the val- 
ues of such that the determinant of the matrix is zero. 
Given the values of E, V, t±, and B this amounts to 
solving a quadratic equation for inside {U (?■) — 0) 
and outside {U{r) = Uq) the QD with the result 



4 




4 



(31) 



which is independent of m. The energies e> and e< are 
defined in Sec. II. With the knowledge of k^, we can 
easily find the corresponding eigenvector '^2- Finally we 
may use Eq. (19) and Eq. (28) to recover the eigenvector 
^' in the original basis. A similar procedure was used 
previously in Ref. 27 in the case of zero magnetic field. 

Given the eigenvectors inside and outside the dot the 
bound state solutions of the full problem are those where 
the two pairs of solutions can be matched at the bound- 
ary of the dot. This is most easily tested by computing 
the determinant of the matrix built up by the four rel- 
evant eigenvectors evaluated at r = i?, where R is the 
radius of the dot. The zeros of the determinant as a 
function of the energy (inside of the gap at r ^ oo) de- 
termines the bound states and their energies. The condi- 
tion of having the determinant equal to zero is the bilayer 
analogue of Eqs. (11) and (12) for the single layer case 
and can straightforwardly be computed numerically al- 
though the analytic expression is long and cumbersome. 



B. Landau levels in biased graphene bilayer 

In this section we briefly review the properties of a 
biased graphene bilayer in a magnetic field. From the 
point of view of the QD there exists one level that is 
of particular importance since it crosses the gap with 
increasing magnetic field (see also Refs. 31,32). In the 
presence of a magnetic field, the Hamiltonian matrix in 
a homogenous system can be written as (r = -1-1) 



fV/2 7at 

7a V/2 

t_L 

\ 



tl 



~V/2 
7a 



7a 

t -V/2/ 



(32) 



in the Landau gauge A = (0, Bx, 0) and for a par- 
ticular sign (s = +1) of the magnetic field (chang- 



spcctrum but with V 
= isy/e\B\/2h[x~ 



—V is obtained). Explicitly, 
i{s/e\B\)(—ihdx)] + ipy with pj, a 



c-numbcr due to translational invariance in y-dircction. 
We have defined 7 = wf^V^/'s- The eigcnstates can 
then be formed by a spinor of the form 



5' = [aAi\n),aBi\n - 1), 0^21"), 052!" + 1)]^ 



(33) 



With this choice the operator matrix in Eq. (32) be- 
comes a matrix of numbers acting on the spinor 4* = 

[aAi,aBi,aA2, aB2]^: 



Ho ~ 



/V/2 7V^ t± 

7V^ V/2 

ti_ - v/2 7V^rTT 

\ 7V^r+T -v/2 J 



(34) 



For n > 1 this leads to a spectrum that (as function of 
7) is very similar to the case without a magnetic field 
as a function of the absolute value of the momentum. 
The most important feature for us is that the gap is still 
present for these quantum numbers. 

For n = —1 the spinor is simply [0, 0, 0, |0)] which leads 
to a flat band (Landau level) at —V/2. 
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FIG. 5: Energy levels in a relatively small bilayer QD (radius 
R = 25 nm) as a function of the magnetic field. The other 
parameters are as follows: t± — 0.4 eV — 15.19hv/R, V — 
1.9hv/R, Uo = 1.52hv/R and s = 1 (i.e. positive B-field). 
The solid and dashed lines are for different valleys. 
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The case n = is more interesting. In this case the 
spinor is of the form [a^ilO), 0, a^2|0), aB2|l)] and the 
resulting problem is the diagonalization of the matrix 

(V/2 \ 

Ho = <± -V/2 7 . (35) 
V 7 'V/2) 

This leads to three levels as a function of 7. Two start 
out at ±t_L at 7 = and evolve smoothly over to ±7 as 
7 ^ 00. These two levels are therefore much like the case 
n > 1. The most interesting level starts out at —V/2 at 
7 = and goes smoothly to 1^/2 as 7 — > cx). It is easy to 
see that the level crosses zero at 7 = \/t\ + V^/A. This 
Landau level that crosses the gap can also be seen in 
Fig. 3 of Ref. 31 and has important consequences for the 
levels in the QD, as wc will discuss in the next subsection. 



C. Results for the bound state levels 

The bilayer QD is in many ways similar to the single 
layer QD discussed above, but there are also important 
differences in the physics. 

The most important result of our study can be seen 
in Fig. 5 where we display the energy levels of a dot 
as a function of the magnetic field. At zero magnetic 
field, the degeneracy of the levels in the two valleys is 
clearly displayed. With increasing the magnetic field, 
the orbital degeneracy is lifted. The symmetry of the 
levels is analogous to the case of the single layer discussed 
above. The states that are degenerate at zero field are 
related by time-reversal symmetry which means that they 
correspond to opposite values of angular momentum ±to 
in different valleys. The typical effective time-reversal 
symmetry of ±m within one valley is already broken by 
a "mass" term (which breaks the inversion symmetry of 
the bilayer) in a similar manner to the case of Neutrino 
billiards considered by Berry and Mondragon.^^ 

An important feature of the bilayer as opposed to the 
single layer is the unconventional Mexican hat-likc dis- 
persion relation near the band edge. This is most appar- 
ent for a large value of the bias field V . An example of the 
level structure for such a dot is shown in Fig. 6. It is clear 
that there are many closely spaced levels near the band 
edge. This is a feature of the enhanced density of states 
near this particular energy. '^'^ It is also crucial to note that 
the trigonal distortion term (which breaks the cylindrical 
symmetry and in principle couples all states with angular 
momenta to, m±3, to±6, to±9, . . .) is a particularly rel- 
evant perturbation for the degenerate states close to the 
band edge. For states away from the band edge for which 
the energies of the coupled states are different in energy, 
the trigonal distortion term can be treated as a pertur- 
bation and we do not expect that the energy levels will 
be much affected. More explicitly, we use a cylindrically 
symmetric dispersion relation whereas the real dispersion 
relation including the trigonal distortion term, which is 
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FIG. 6: Energy levels in a relatively small bilayer QD with 
the same parameters as in Fig. 5 except that the bias field V 
is about three times larger; V — 6hv/R. 



parametrized by V3 (v^ w 0.1 in graphite but it has not 
yet been measured in bilayer graphcne), is not. From 
the expressions in Ref. 27 (zero magnetic field) we find 
that above the momentum scale Pc ~ V3t± the cylindri- 
cally symmetric term that we keep is dominating over 
the trigonal distortion term in the Hamiltonian and does 
hence provide a reasonable zeroth order approximation. 
It is not trivial to convert this momentum scale into an 
energy in general because of the Mexican hat structure. 
But for the parameters we use in Figures 5 , 6 and 7 the 
associated momentum is larger than for energies above 
the region of the Mexican hat (i.e. above V/2) where the 
energy becomes a monotonously growing function of mo- 
mentum. Therefore, we believe that our model captures 
the relevant physics above the Mexican hat. We note 
that at finite magnetic fields, it is known that the trigo- 
nal distortion quickly becomes less important with grow- 
ing magnetic field in an unbiased bilayer. '^^ Wc therefore 
expect that for the large field regime, the corrections are 
small at all energies. Additional subleading parameters 
(such as 74 which introduce an electron-hole asymme- 
try into the spectrum^^) will also shift the level positions 
slightly. 

For a large QD, it is also possible to reach the regime 
where the dot levels are described by the Landau lev- 
els. This feature is seen in Fig. 7a) where wc display 
the bound states for to = 0,to = ±1 for large magnetic 
fields. Note that the QD levels tend to approach the bulk 
Landau levels displayed in Fig. 7b). In a smaller QD it is 
hard to reach the Landau level limit for moderate mag- 
netic fields. 

Another important feature of the bilayer for designing 
a QD is the existence of an anomalous LL that crosses 
the gap, see subsection III.B and Fig. 7b). The character 
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FIG. 7: a): Merging of bilayer QD levels to the bulk LLs as 
function of magnetic field for a relatively large bilayer QD 
with R=67.48 nm and t± = 0.4 eV = Alhv/R, Uo = 3.5Hv/R, 
V = 5.13hv/R for m = 0, ±1 and s = 1 (i.e. positive B- 
field). Full lines are for r = +1 and dashed lines are for 
r = — 1. b) Bulk Landau Levels (LL) which are approached 
almost perfectly at high fields in this parameter regime. Note 
that the n = LL (see subsection III B) crosses the gap with 
a negative slope, whereas the other LLs (n — 1,2,3) have 
positive slope (blue and red denotes the two valleys). There 
is also a flat LL (n = —1) at V/2 similar to the single layer 
case. 

of bound states changes, when the square root of ^ in 
Eq. (31) changes sign, which occurs at energies 

srV I VHl t\ 

< i\{t\^v^) ^m+v) l%itl+V^r' 

(36) 

and at Ey = i?< + Uq. These lines are shown in Fig. 8 
as a function of magnetic field. The area between i?< 
and define an effective (valley-dependent) bandwidth 
for the QD. Indeed, at _B = 0, we obtain only bound 



states for energies above \t±V\/2y/t'^ + and below 
\t±V\/2y^^^+V^ +Uo which correspond to the conduc- 
tion band minimas inside and outside the QD, respec- 
tively. Within this bandwidth, the two k^'s inside the 
QD (k<) are purely real and the two k^'s outside the 
QD (k>) have an imaginary part. The physical meaning 
of n/ {\/21b) is most transparent in the limit of zero mag- 
netic field where it becomes the inverse decay length of 
the wave function. ■^^ Thus states within the bandwidth 
correspond to a decaying wave outside that is matched 
to one propagating and one decaying wave inside of the 
dot (this is true for energies such that only one band is 
allowed inside of the dot whereas the other band requires 
the momentum to be imaginary). 

The bound state energy window at B = is crossed by 
the n = bulk LL, as shown in Fig. 8. When QD bound 
states cross this bulk LL, the QD becomes "leaky" and 
electrons can escape into the bulk. The effective band- 
width defined by the area between the lines -E> and i?< is 
never crossed by a bulk LL and defines therefore a " safe" 
zone for QD bound states. Note, however, that bound 
states do exist also outside of this effective bandwidth at 
finite magnetic field, since the the bulk spectrum (LLs) 
becomes discrete. This means that evanescent waves con- 
tinue to exist in the " bulk" region when they are not de- 
generate with a bulk LL (much like the edge states that 
are present between the LLs in the integer Quantum Hall 
Effect). However, this regime is not ideal for QDs, due 
to the leakage via nearby bulk states as described above. 

We point out that QDs in bilayer graphene in con- 
nection with a magnetic field again allow for a controlled 
tuning of level degeneracies. The values used in our plots 
correspond to realistic values of the gap voltage V and the 
inter-layer coupling t±}^ Besides similarities to the sin- 
gle layer case studied in section II, the bilayer QD shows 
very interesting additional features. The size of the gap 
V can influence the size of the level spacing drastically 
when the energies are close to the band edge, where a 
Mexican-hat like structure is formed which is more pro- 
nounced at larger V . 

IV. APPLICATIONS OF VALLEY SPLITTINGS 

In this section we discuss the implications and their 
use of the broken valley degeneracy by a magnetic field 
in gate-tunable graphene QDs. 

A. Consequences for valleytronics 

The valley index r = ±1 can be thought of as eigenval- 
ues of the operator u ■ r where i/ is a unit vector on the 
Bloch sphere and r the vector of Pauli matrices. The 
operator r is called the valley isospin. If the two val- 
leys arc uncoupled, we have u = z. It has been pointed 
out that the valley isospin r could be used as a con- 
trollable degree of freedom like the electron spin S is 
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FIG. 8: Effective bandwidth of the bilayer QD defined by the 
area between the lines iJ> and (dashed lines) as a function 
of magnetic field for V = VJ.lhv/R and Ua = 7hv/R, s = 1, 
t± — Alhv/R and for valley t = — 1. The dotted line displays 
the n = bulk LL (in valley r = — 1) that crosses the QD 
bandwidth corresponding to the band gap at zero magnetic 
field (full lines). This bulk LL presents an escape channel for 
the bilayer QD if bound states cross this bulk LL. Note that 
no bulk states overlap with the effective bandwidth (from the 
same valley). 



used in spintronics applications which coined the name 
valleytronics.^^ The main motivation to use the valley de- 
gree of freedom as a new unit of information in graphene 
instead of the sublattice pseudospin cr, is the fact that 
the valley degree of freedom is preserved in the absence 
of short range scatterers (whereas cr is not), provided 
e.g. by the graphene edges. However, the manipulation 
of the valley isospin is not as straightforward as for the 
real electron spin since the valley isospin does not di- 
rectly couple to a magnetic field as does the real spin via 
the Zeeman interaction. However, since the valleys are 
related by time-reversal symmetry, the valley degeneracy 
can also be broken in principle by applying a magnetic 
field which we have shown in this work. However, a mag- 
netic field alone is not enough since it breaks only degen- 
eracies within different valleys. The so called effective 
time reversal symmetry p —f —p and cr — s- — cr within 
each valley should also be broken. This is achieved by 
quantum confinement induced by a boundary that does 
not couple the valleys,^" which is the case for the gate- 
tunable QDs proposed here (see also Sec. II B). 

In this work we have shown that the valley degeneracy, 
and more generally, the orbital degeneracy is controUably 
and efficiently broken by a magnetic field. In the case of 
the valley splitting Ak,k', we can take advantage of the 
anomalous LLs which are approached by the QD states 
at large R/Ib- As we have discussed in sections II and 
HI, there exists a flat LL at the gap value that is only 



present in one of the valleys (its other valley partner is at 
negative energy, i.e. in the valence band). In the bilayer 
QD, we have in addition a state that crosses the gap with 
increasing magnetic field. We can estimate typical values 
for l^K,K' by comparison with our plots for QD bound 
states. In the single layer, we obtain from Fig. 3b) for a 
dot radius of i? = 67.48 nm and for B ~ 4.6 T, a valley 
splitting /S.K.K' of about 24.4 meV between the valley- 
polarized ground state and the first excited states (from 
the other valley). For the same QD size and for the same 
S-field, we obtain for the bilayer QD shown in Fig. 7a), 
a splitting between the approached n = LL and the 
n = 1 LL from the other valley of Ak,k ~ 8 meV. Both 
values for Ak,k' are much larger than temperature and 
can be probed in timneling transport through QDs. 

Such QDs could be used as very efficient valley-filters 
in transport through the QD. In contrast to earlier pro- 
posals of valley filters in zigzag ribbons in single layer 
graphene, and topologically confined states in bilayer 
graphene'^^ without magnetic field, the present setups 
would work as filters that break time-reversal symmetry 
and therefore function also in closed QD systems where 
Coulomb blockade (CB) effects can be used to operate at 
the single (valley-) spin level much like for ordinary spin- 
filters in QDs."^^ CB effects (and therefore single electron 
tunneling) become prominent if the charging energy ex- 
ceeds the temperature and for weak coupling to the leads 
(tunnel resistance ^ h/e^). We estimate the charging en- 
ergy Eq as a function of the QD radius R as Eq = e^/C 
with capacitance C = SeoScsR where EcS = (1 + 4)/2, 
including the dielectric constant for Si02 and vacuum.® 
For R = 67.48 nm, this gives Ec ^ 12 meV. The tunnel- 
ing rate in and out of the QDs could be tuned by gates. 
We note that two valley-filters in series can be used as a 
valley valve,^^ where the valley polarization of one of the 
QDs should be reversible. This can be achieved by either 
reversing the sign of the magnetic field, or more easily, by 
gates such that one QD can be tuned from a hole-doped 
QD to a n-doped QD (and vice versa). In this way, the 
valley isospin of QD states at resonance with the the 
leads can be changed in one QD, thereby probing the po- 
larization of the other QD. We note that the presence of 
such a valley splitting could be probed by electron trans- 
port since the level degeneracy is changed with increasing 
magnetic field. Since the valley splitting Ak,k' acts like 
a Zeeman field for the valley isospin, experiments that 
measured the relaxation time Ti ,'^® and the read-out of a 
single electron spin'^^ in GaAs QDs could be performed 
in a similar way in gate-tunable graphene QDs in order 
to measure the valley relaxation time and valley polar- 
ization in such QDs (besides the detection measurements 
for real spin). 



B. Consequences for spin qubits 

The use of the spin 1/2 degree of freedom of single 
electrons as qubits^" is usually combined with a proposed 
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coupling of adjacent localized spin qubits via the Heiscn- 
berg exchange interaction.*^ Carbon may offer relatively 
long spin coherence times due to the sparseness of nu- 
clear spins and potentially also due to the weakness of 
the spin-orbit coupling in these materials. However, 
spin qubits in graphenc QDs,^^ need to deal with the 
valley degeneracy in these materials which can interfere 
with exchange coupling. This can be understood using 
the following simple model for two electrons in adjacent 
graphene QDs. Suppose that the valley degree of free- 
dom, unlike the spin in this case, is not well under con- 
trol, and, because it is degenerate, each electron is in an 
incoherent mixture of the two valley states K and K', 
with equal probability, pr = {\K){K\ + \K'){K'\)/2. The 
density matrix of the two electrons in adjacent QDs can 
thus be written as pr = {\KK){KK\ + \K'K){K'K\ + 
\KK'){KK'\ + \K'K'){K'K'\)/4. Including spin, the 
density matrix is then p = pr ® \(p) {<p\ where \ip) is an 
arbitrary pure two-spin state. At this point, the spin and 
valley degrees of freedom are uncoupled, and while the 
spin can maintain its coherence, the valley isospin may 
at the same time be in an entirely incoherent state. The 
problem arises because there will be a tunnel-coupling 
mediated exchange coupling J 7^ if the electrons arc 
both in the same valley \KK) or \K'K') but there will 
be no such coupling (J = 0) in the cases where the elec- 
trons are in different valleys, i.e. \KK') and \K'K). The 
reason for this is that the exchange coupling relies on the 
Pauli exclusion principle which only matters in case that 
both electrons can occupy the same orbital. Here, we 
assume that the inter-dot tunneling conserves the val- 
ley isospin (however, we note that a similar conclusion 
would be obtained in the non-conserving case). Suppose 
we apply the exchange coupling such that it generates 
a SWAP operation that exchanges the two spin qubit 
states. *° This SWAP operation will be conditional on the 
valley state. Assuming an initial spin state = | H — ), 
where |±) = (| t) ± | i))/V2, we find, after the SWAP, 
the state p' = {\KK){KK\ + \K'K'){K'K'\)/2(x)\- +) + 
{\KK'){KK'\ + \K'K){K'K\)/2 ® | + -). If the valley 
degree is traced out, we find that p' ~ {\ — h)( — \- \ + 
I H — )(H — |)/2. With this, the phase coherence of both 
spins decays (cr^) = 0, i = 1, 2, due to the coupling to the 
incoherent valley degree of freedom. Even if the valley 
isospin is coherent, a valley degeneracy will still lead to 
spin-isospin entanglement, which for some purposes may 
be interesting, but which essentially reduces the spin co- 



herence to some charge (valley) coherence time which 
can be expected to be shorter. If the valley degeneracy 
can be lifted, as proposed here, one can avoid this en- 
tanglement and possible spin decoherence processes that 
arise from it. We note that orbital degeneracies in the 
same valley lead to similar problems for the exchange 
coupling of neighboring spins. We therefore should oper- 
ate at moderate magnetic fields such that the LL regime 
is not reached. 

V. CONCLUSIONS 

We have studied the bound states of QDs in gapped 
single- and bilayer graphene in the presence of a homo- 
geneous magnetic field. Due to the absence of sharp 
graphene edges, the valleys are well defined in these QDs. 
We have shown that these realistic structures would al- 
low us to control the valley degeneracy by the magnetic 
field. This has important consequences for spin or valley- 
quantum computing, where the breaking of orbital (or 
valley) degeneracy is absolutely crucial. Besides similar- 
ities between the two systems, we also found crucial dif- 
ferences that can be traced back to an anomalous Landau 
level (LL) in the gapped bilayer that crosses the gap and 
which can provide an escape channel for QD bound states 
into the bulk at large magnetic fields if they cross this LL. 
In addition, the level spacing size close to the band edge, 
crucially depends on the strength of the applied voltage 
in the bilayer QD, which is due to a "Mexican hat" form 
of the bulk bandstructure. These features have impor- 
tant implications for finding the ideal parameter range 
for useful QDs. We also discussed possible applications 
of such QDs for the emerging fields of vallcytronics and 
spin qubits in graphene. 
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